#include <stdio.h>
#include <time.h>
#include <limits.h>
#include <gmp.h> // gcc -lgmp
#include <pthread.h>

#define N 25
#define BASE 10
#define MANTISSA 1024*1024*32

#define N_THREADS 2+(N*4)

/*
N = 10.000
MANTISSA de 128K
104.00 segundos

*/

void *a_zero(); // OK = 0.34314...
void *y_zero(); // OK = 0.41421...
void *eq_exp(void *n);
void *eq_ay4();
void *eq_yep();
void *y_kn(); // Ok 0.3734i...^e-2

mpf_t expt;
char isExpInit; // Serve para verificar se exp ja foi inicializado

mpf_t ay4;
char isAy4Init; // Serve para verificar se ay4 ja foi inicializado

mpf_t yep;
char isYepInit; // Serve para verificar se yep ja foi inicializado

mpf_t a;
mpf_t y;

int main(int argc, char* argv[]){

	unsigned long int k=0, i=0, t=0;
	time_t t_zero, t_final;
	double delta_t;
    
	mpf_t pi;

    pthread_t threads[N_THREADS];
    pthread_attr_t attr;

	time(&t_zero);

    pthread_attr_init(&attr);
    pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);

	mpf_set_default_prec(MANTISSA);

    isExpInit=0;
    isAy4Init=0;
    isYepInit=0;
	
    // mpf_t a && mpf_t y são inicializdos dentro das threads
    pthread_create(&threads[i], &attr, a_zero, NULL); i++;
    pthread_create(&threads[i], &attr, y_zero, NULL); i++;

    pthread_join(threads[t], NULL); t++;
    pthread_join(threads[t], NULL); t++;
    
	/*

	a = 6 - (4 * 2^1/2)
	y = (2^1/2) - 1

	yn = (1-(1-y^4)^1/4)/
	     (1+(1-y^4)^1/4)

	an = a(1+ yn)^4 - 2^(2k+3) ? yn(1 + yn + yn^2)
	**	an converge para 1/pi

	*/
	
	for(k=0; k<N; k++){// K na primeira interação é = 0

        pthread_create(&threads[i], &attr, y_kn, NULL); i++;
        pthread_create(&threads[i], &attr, eq_exp, (void*)k); i++;

        pthread_join(threads[t], NULL); t++; // espero a y_kn
        
        pthread_create(&threads[i], &attr, eq_ay4, NULL); i++;
 
        pthread_join(threads[t], NULL); t++; // espero a eq_exp
 
        pthread_create(&threads[i], &attr, eq_yep, NULL); i++;

        pthread_join(threads[t], NULL); t++; // espero eq_ay4
        pthread_join(threads[t], NULL); t++; // espero a eq_yep

        mpf_sub(a, ay4, yep);    
    }

	mpf_clear(y);
	mpf_clear(yep);
	mpf_clear(ay4);
	mpf_clear(expt);

    pthread_attr_destroy(&attr);

	mpf_ui_div(a, 1, a); // an converge para 1/pi
	mpf_out_str(NULL, BASE, 0, a);
	mpf_clear(a);
	
	time(&t_final);
	delta_t = difftime(t_final, t_zero);
	printf("\nTempo = %.2lf segundos\n", delta_t);	

    pthread_exit(NULL);

	return 0;
}

void *a_zero(){ // OK = 0.34314...

	mpf_init(a);
	// a = 6.0-(4.0*sqrt(2));
	mpf_sqrt_ui(a, 2); 
	mpf_mul_ui(a, a, 4);
	mpf_ui_sub(a, 6, a);

    pthread_exit(NULL);
}

void *y_zero(){ // OK = 0.41421...

	mpf_init(y);
	//y = sqrt(2)-1.0;	
	mpf_sqrt_ui(y, 2); 
	mpf_sub_ui(y, y, 1);

    pthread_exit(NULL);
}

void *eq_exp(void *n){
    
    unsigned long int k;

    if(!isExpInit){
        mpf_init(expt);
        isExpInit=1;
    }

    k = (unsigned long int) n;
     
    k= (2*k)+3;
	mpf_set_ui(expt, 2);
	mpf_pow_ui(expt, expt, k);	

    pthread_exit(NULL);
}

void *y_kn(){ // Ok 0.3734i...^e-2

	mpf_t y4, d;	
	mpf_init(y4);	
	mpf_init(d);
	
	//	yn = (1-(1-y^4)^1/4)/
	//	     (1+(1-y^4)^1/4)

	mpf_pow_ui(y4, y, 4);
	mpf_ui_sub(y4, 1, y4);
	mpf_sqrt(y4, y4);
	mpf_sqrt(y4, y4);
	
	mpf_add_ui(d, y4, 1);
	mpf_ui_sub(y4, 1, y4);

	mpf_div(y, y4, d);

	mpf_clear(y4);
	mpf_clear(d);

    pthread_exit(NULL);
}

void *eq_ay4(){
    
    if(!isAy4Init){
        mpf_init(ay4);
        isAy4Init=1;
    }
 
	mpf_add_ui(ay4, y, 1);
	mpf_pow_ui(ay4, ay4, 4);
	mpf_mul(ay4, a, ay4);

    pthread_exit(NULL);
}

void *eq_yep(){

    mpf_t yk;

    if(!isYepInit){
        mpf_init(yep);
        isYepInit=1;
    }
 
	// yk = (1+yn+yn^2)
	mpf_init(yk);
	mpf_pow_ui(yk, y, 2);
	mpf_add(yk, yk, y);
	mpf_add_ui(yk, yk, 1);

    // yk = yn(yk)
	mpf_mul(yk, yk, y);	

	// exp*yk
	mpf_mul(yep, expt, yk);

    mpf_clear(yk);

    pthread_exit(NULL);
}
